The ferromagnetic transition and domain structure in LiHoF4 
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Using Monte Carlo simulations we confirm that the rare-earth compound LiHoF4 is a very good 
realization of a dipolar Ising model. With only one free parameter our calculations for the mag- 
netization, specific heat and inverse susceptibility match experimental data at a quantitative level 
in the single Kelvin temperature range, including the ferromagnetic transition at 1.53 K. Using 
parallel tempering methods and reaching system sizes up to 32000 dipoles with periodic boundary 
conditions we are able to give strong direct evidence of the logarithmic corrections predicted in 
renormalization group theory. Due to the long range and angular dependence of the dipolar model 
sample shape and domains play a crucial role in the ordered state. We go beyond Griffiths's theorem 
and consider surface corrections arising in finite macroscopic samples leading to a theory of magnetic 
domains. We predict that the ground-state domain structure for cylinders with a demagnetization 
factor A*' > consists of thin parallel sheets of opposite magnetization, with a width depending on 
the demagnetization factor. 

PACS numbers: 75.10.Hk,75.40.Cx,75.40.Mg,75.50.Dd,75.60.Ch 



The use of effective theories is one of the primary 
modus operandi of modern physics. Frequently the effec- 
tive models give only a qualitatively accurate description 
of the phenomena under investigation, due to corrections 
that are omitted and free parameters that may be hard 
to determine experimentally. Finding experimental sys- 
tems that display striking phenomena, are accurately de- 
scribed by a simple model, and have few or no free param- 
eters, is important since it enables detailed comparison 
between experiments, theory and numerical simulations. 

The rare earth magnet LiHoj,Yi_a:F4 displays an ar- 
ray of fascinating magnetic phenomena such as quantum 
phase transitionSfi spin-glass behavior,^ and persistent 
coherent oscillations^. Yet at least the pure material 
LiHoF4 is believed to be described by one of the most 
fundamental models in condensed matter physics: the 
two-state Ising model. Materials such as the antifer- 
romagnets DyP04 and DyAlG have been shown to be 
accurately described by a short-ranged Ising model4. In 
LiHoF4, on the other hand, the magnetic properties are 
dominated by the dipolar interaction. Since the inter- 
action strength is set by the known g factor this makes 
it possible to determine the effective model to high accu- 
racy. However, the inherent frustration and long-range of 
the dipolar model make direct numerical simulations de- 
manding. Using a parallel tempering Monte Carlo (MC) 
method that is essentially free of systematic errors (apart 
from finite-size effects) we go beyond mean-field theory 
and explicitly demonstrate that the experimental data 
for LiHoF4 is indeed in quantitative agreement with the 
dipolar Ising model. 

The magnetic properties of LiHoF4 originate in the 
4/-electrons of the Ho^+ ions which sit in a tetrago- 
nal lattice with a unit cell of size (1,1,2.077) in units 
of a = 5.175^. According to Hund's rules the holmium 
ion has a ^1% ground state, but the crystal field partially 
lifts the 17-fold degeneracy, and the resulting doubly- 
degenerate ground state is separated from the first ex- 



cited state by 11 K. This separation of energy levels en- 
ables a projection of the full Hamiltonian onto the ground 
state subspace.^ The matrix elements of the operators 
and jy vanish in this subspace, and the effective model 
is the dipolar Ising model 

The dipolar coupling constant is given by Jd = 
{qI^b /"^Y / 0? — 0.214 K due to the renormalized g fac- 
tor ~ 13.8, which can be computed from the crystal-field 
Hamiltonian^, or deduced from the experimental high- 
temperature susceptibility^. The only free parameter in 
the model, the weak exchange interaction, has no funda- 
mental physical effect on the system other than to alter 
the critical temperature (Tc)^ii. We set it to Je = 0.12K 
to reduce the Tc of the model from 1.91 K to the experi- 
mental value 1.53 K. 

The study of the dipolar interaction has a long and 
interesting history. Due to demagnetization effects, Lut- 
tinger and Tisza* found a ground state energy that de- 
pends on both lattice structure and sample shape. Grif- 
fiths later gave a proof that the free energy, without ap- 
plied fields, is independent of sample shape. ^ The appar- 
ent contradiction is explained in terms of domain forma- 
tion, allowing the magnetic order to vary from one macro- 
scopic part of the system to the next. Experimentally this 
has been demonstrated since measurements of the specific 
heat for LiHoF4 show no apparent shape dependence^ 
and needle shaped domains have been observed close to 
the transition. ^-"^ 

The dipolar interaction has several properties that 
complicate a numerical treatment of the model. Inher- 
ent frustration combined with the long range makes MC 
equilibration cumbersome at low temperatures, requir- 
ing long simulation runs to reach equilibrium. To handle 
the long range of the interactions we employ the method 
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of Ewald summationi^.. This method not only gives im- 
proved numerical convergence due to the use of periodic 
boundary conditions, but also includes an Ewald param- 
eter which emulates different sample shapes. Our sim- 
ulations have been carried out using single-flip parallel 
tempering MC, since cluster methods are of limited use 
in frustrated systems. The MC sample is of size Li^ unit 
cells, with the linear size L ranging from 10 to 20 (4000 
to 32000 spins) throughout the study. In most figures we 
use about 100 temperature points resulting in smooth 
curves. 




FIG. 1: The internal energy for spherical (black lines) and 
cylindrical (dashed red line) boundary conditions. For the 
spherical boundary the MC cell is of size unit cells with 
L = 10, 12, 14, 16 and 18 from top to bottom. The inset shows 
the difference of the two energies as a function of system size 
at T=0.5K. 



In order to investigate the convergence to the ther- 
modynamic limit we first consider the effects of different 
sample shapes on the internal energy in Fig. ([1]). The 
calculation is performed for a long needle (with demag- 
netization factor N=0) and for a sphere (N=47r/3). The 
energy for the needle converges quickly in system size and 
we show the converged curve. The energy for the sphere 
coincides with the energy for the needle above Tc, but 
shows large finite-size corrections below Tc- The needle 
orders ferromagnetically and according to Griffiths 's the- 
orem the infinite spherical sample must have the same 
energy. If the spherical sample forms ferromagnetic do- 
mains that cancel the internal magnetic field, then the 
two energies will be equal. The formation of domains 
can be seen directly in the simulations of the spherical 
sample, and in Fig. ^ we show a typical spin configura- 
tion at T = O.STc. In the inset of Fig. ^ we see that the 
difference of the two energies decreases as the inverse of 
the volume of the system. 

In order to verify the accuracy of our effective model 
for LiIIoF4 we make detailed comparisons between our 
calculations and existing experimental data. In Fig. (|3l) 
we show the specific heat measurements from Refs. llOlflil 
and our calculations for spherical and cylindrical (N=0) 
samples. Again we find that the numerical results for 





FIG. 2: Typical spin configuration for spherical boundary 
conditions at T = O.STc. 



MC cylinder 

MC sphere 

♦ fromRef. [10] 
■ fromRef. [13] 




FIG. 3: The specific heat capacity for spherical (L=12,10 and 
8 from top to bottom) and cylindrical boundary conditions 
(L=20 and 18 from top to bottom), and experimental data 
for a spherical samplei^ as well as an oblate sample^*^. 



spherical samples show slow convergence, while the re- 
sults for zero demagnetization factor have converged, ex- 
cept for close to the critical temperature, where finite-size 
effects are still visible. In the range of converged data the 
numerical results agree very well with the experimental 
data. 

Next we compare the spontaneous domain magnetiza- 
tion measured by Griffin et al.^^ with our simulations for 
zero demagnetization factor in Fig. The experimen- 
tal data is only determined up to a constant, and we have 
normalized the experimental data to agree with our cal- 
culations at 1.35 K. The agreement is very good all the 
way up to about 0.96Tc, where finite size effects become 
visible in our largest system sizes. Below we will analyze 
our data more carefully and demonstrate that we can di- 
rectly observe logarithmic corrections to the mean-field 
magnetization. 

Finally we consider the inverse susceptibility as mea- 
sured by Cooke et al^ as a function of different de- 
magnetization factors. In Fig. ^ we demonstrate that 
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FIG. 4: Magnetization as a function of temperature for sys- 
tem sizes L=10,18 and 20, as well as experimental data^'*. 
Close to the critical temperature we show a fit to our data 
including logarithmic corrections. 



the form log|(T - T^) /Tq\'^I'^ , where Tq is an effective 
temperatureii^ Experimentally, logarithmic corrections 
have been convincingly seen in the magnetizationi^ of 
LiHoF4 and the specific heat of LiTbF4^^, but not in the 
susceptibility.^^ Numerical studies of the dipolar model 
have applied finite-size scaling to detect the logarithmic 
corrections, ^^'^^ but since we have results for large system 
sizes (32000 spins) we directly fit a curve of the form 

m{t) ^ (T - T,)i/2| log |(r „ Te)/To|r (2) 

to the part of the critical region where the MC data has 
converged. For different values of Tc we let the Tq and 
the exponent a vary and display the value of a that gives 
the best fit, together with the corresponding value 
in Fig. ([5]). There is a minimum in around a = 0.18, 
giving numerical evidence of non-zero logarithmic correc- 
tions. However, the exact value of the optimal and a 
depends on the temperature interval included in the fit, 
but a finite value of a does improve the fit. 



the agreement with experiments again is very good for 
N = 47r/3 (spherical sample) and N=1.65. The inter- 
esting behavior of the the inverse susceptibility below Tc 
can be understood in terms of the domains ;S With no 
disorder the walls are free to move and arrange them- 
selves to cancel the internal field: = H — A^M = 0, 
resulting in a constant susceptibility, x = M/H = 
below Tc. We find it remarkable that the numerical sim- 
ulations for limited system sizes have converged so well, 
even if the domain size is much smaller than in the real 
material. Above Tc the Curie- Weiss law is followed. For 
the N = 1.65 sample we note that the initial slope of the 
experimental data is slightly higher than the MC data, 
indicative of a higher g factor. 
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FIG. 5: The inverse magnetic susceptibility 1/x from experi- 
ments and simulations. The three sets of curves correspond to 
N=47r/3 (spherical sample), 1.65 and from top to bottom. 

The upper critical dimension is three for uniaxial 
dipolar interactions, and according to renormalization 
group theory, the magnetization, susceptibility and spe- 
cific heat are predicted to have logarithmic corrections of 
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FIG. 6: Exponent a of the logarithmic correction and x for 
optimal fit of the magnetization curve as a function of Tc. 

For the thermodynamic quantities that diverge at the 
critical point we have not been able to use the same direct 
fit due to the large finite-size effects. However, convincing 
evidence for logarithmic corrections in the heat capacity 
can still be obtained by plotting the peak height against a 
logarithmically corrected system size. This curve should 
tend to a constant for large system sizes, and as can be 
seen in Fig. ([7]) the curve levels out significantly faster for 
the predicted exponent a 1/3 than for the mean-field 
result a — Q. 

Griffiths's theorem predicts the formation of domains 
in order to make the free energy independent of sample 
shape, but it does not answer the fundamental question of 
the size and shape of the domains. The domain structure 
is the result of an energy balance: introducing a domain 
wall increases the dipolar energy of Eq. ([T|) , whereas the 
magnetostatic energy^ is decreased. The energy per spin 
is of the form E = Ci/ D + C2 ■ D + C-i, where D is the 
linear size of the domain, Ci the domain wall energy, 
C2 the magnetostatic energy and C3 is the energy of the 
ferromagnetic ground state for = 0. 
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FIG. 7: The peak of the heat capacity grows logarithmically 
with system size. 





FIG. 8: Domain configurations of parallel sheets (-E2) 
checkerboard pattern (-E4) and cylinders (Ec). 



The magnetostatic energy arises from the surface 
charge of the domain configuration. In the limit of in- 
finite system size this term vanishes, and the energy is 
given by E = Ci/ D + C3, which is the hmit considered 
by Griffiths. We consequently see that as the domain 
size increases the energy approaches C3. The energy of 
a particular domain configuration can be calculated as a 
function of D using the Ewald summation, enabling us 
to extract the constant Ci for any domain configuration. 

Here we consider three domain structures that are 
likely candidates for the ground state of LiHoF4: thin 
parallel sheets {E2), checker board {E4) and cylinders 
{Ec), depicted in Fig. ([5]). Kittel has calculated the cor- 
responding magnetostatic surface energy density as 0.85, 
0.53 and 0.37 ^APD?° In Fig. (O we plot the energy per 
spin as a function of domain size in the fully polarized 
ground state for the three configurations. The configu- 
ration of parallel sheets has the lowest energy, and we 
predict that this is the ground state domain structure. 
The calculation was done for a cylinder of diameter 3.2 



mm and length 4.8 mm. We also show E2 for lengths 
9.6 mm and 19.2 mm. As the demagnetization factor 
of the cylinder decreases, the size of the domains grows, 
as shown in Fig. Q, but E2 remains the lowest energy. 
From Fig. ([9]) we see that for finite samples there is ac- 
tually a small shape dependence of the energy, which 
disappears in the limit of infinite system size considered 
by Grifliths. 

We have provided strong evidence that the rare-earth 
magnet LiIIoF4 is a very good realization of the dipo- 
lar Ising model. This enables detailed comparisons of 
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FIG. 9: The energy per spin for parallel sheets (-E2), checker- 
board (-E4) and cylindrical (Ec) domain structures for a cylin- 
der of diameter 3.2 mm and height 4.8 mm (solid curves). 
Dashed curves show E2 for cylinders of increasing heights 9.6 
mm and 19.2 mm. 



theory, experiments and simulations. As examples of 
this we verify the logarithmic corrections predicted by 
renormalization group theory for the dipolar model, and 
we predict that the ground state domain configuration 
for cylindrically shaped samples consists of thin parallel 
sheets. Since the domain walls have no width, very clean 
single crystals with no signs of strain are available, and 
domains appear naturally in MC simulations, LiIIoF4 is 
an excellent testing ground for theories of domains. In 
particular we believe there is much further scope for the 
study of domain-wall motion in the presence of disorder 
and transverse fields. 
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^ D. Bitko, T. F. Rosenbaum, and G. Aeppli, Phys. Rev. 

Lett. 77, 940 (1996). 
^ D. H. Reich, B. EUman, J. Yang, T. F. Rosenbaum, 

G. Aepph, and D. P. Belanger, Phys. Rev. B 42, 4631 

(1990). 

^ S. Ghosh, R. Parthasarathy, T. F. Rosenbaum, and 
G. Aepph, Science 296, 2195 (2002). 



W. P. Wolf, Braz. Jour. Phys. 30, 794 (2000). 
^ P. B. Chakraborty, P. Henelius, H. Kj0nsberg, A. W. Sand- 

vik, and S. M. Girvin, Phys. Rev. B 70, 144411 (2004). 
^ A. H. Cooke, D. A. Jones, J. F. A. Silva, and M. R. Wells, 

J. Phys. C 8, 4083 (1975). 

A. Bihmo and P. Henehus, Phys. Rev. B 76, 054423 
(2007). 



5 



J. M. Luttingcr and L. Tisza, Phys. Rev. 70, 954 (1946). 
R. Griffiths, Phys. Rev 176, 655 (1968). 
G. Mcnncnga, L. ,J. dc Jongh, and W. J. Huiskamp, J. 
Magn. Magn. Mater. 44, 59 (1984). 

J. E. Battison, A. Kasten, M. J. M. Leask, J. B. Lowry, 

and B. M. Wanklyn, J. Phys. C 8, 4089 (1975). 

P. P. Ewald, Annalen der Physik 369, 253 (1921). 

J. Nikkei and B. Ellrnan, Phys. Rev. B 64, 214420 (2001). 

J. A. Griffin, M. Huster, and R. J. Folweiler, Phys. Rev. B 

22, 4370 (1980). 

A. I. Larkin and D. E. Khmel'mtskij, Soviet Physics JETP 



29, 1123 (1969). 

G. Ahlers, A. Kornblit, and H. Guggenheim, Phys. Rev. 
Lett. 34, 1227 (1975). 

P. Beauvillain, J. P. Renard, I. Laursen, and P. J. Walker, 
Phys. Rev. B 18, 3360 (1978). 

H. -J. Xu, B. Bergersen, and Z. Racz, J. Phys.: Cond. Mat. 
4, 2035 (1992). 

A. V. Klopper, U. K. Rofiler, and R. L. Stamps, Eur. Phys. 
J. B 50, 45 (2006). 

C. Kittel, Rev. Mod. Phys. 21, 541 (1949). 



